Introduction

Dynamical systems are often defined in terms of the variation we’re experiencing given the point we are in:

\[ \dot{X} = f(X) \]

An increase in the number of the system components, as well as in the interaction complexity, easily makes the system chaotic, showing unpredictability in the long run and extreme sensitivity on the initial conditions. If we keep increasing the system complexity it will eventually become fully random, for instance this is what happens in the case of particles motion whose position can be thought as the sum of iid random variables. This kind of motion is named Brownian Motion.

Our project aims to investigate the price evolution in the stock market eventually proving if it can be considered a form of random walk, specifically a Geometric Brownian Motion (GBM), unpredictable by definition, or if it can be thought as a Chaotic System instead.

To do so we’ll use various methods used to characterize dynamical systems specifically:

We’ll then compare the results coming from a simulated system with the ones from real world data and we will conclude reporting if any difference can be appreciated.

Random Walks - Geometric Brownian Motion

Before proceeding we first have to define what Random Walks and Geometric Brownian Motion are.

Random Walks

A random walk is a random process that describes a path that consists of a succession of random steps. Below we can see an example of 2D random walk on \(\mathbb{Z}^2\) which start at \(0,0\) and at each time step the point moves \(-1\) or \(1\) with equal probability on each axis.

Geometric Brownian Motion

The Brownian motion is the random motion of particles suspended in a medium, it can be described mathematically by the Wiener process \(W_t\), a continuous-time stochastic process characterized by the following properties:

  • \(W_0 = 0\);
  • \(W_t\) is almost surely continuous;
  • \(W_t\) has indepedent increments;
  • \(W_t - W_s \sim \mathcal{N}(0,\ t-s)\) for \(0\leq s \leq t\), where \(\mathcal{N}(\mu,\ \sigma)\) is the normal distribution with expected value \(\mu\) and variance \(\sigma^2\).

The Wiener process can be constructed as limit of a random walk.

A geometric Brownian motion (GBM) is a continuous-time stochastic process in which the logarithm of the process follows a Brownian motion, often with a drift. It can be defined as the process that satisfies the following stochastic differential equation (SDE):

\[ d S_t = \mu S_t d t + \sigma S_t dW_t \] where \(W_t\) is a Wiener process, \(\mu\) is the drift and \(\sigma\) is the amount of volatility.

Rewriting the former SDE into

\[ \frac{d S_t}{S_t} = \mu dt + \sigma dW_t \]

we can interpret it as the percentage change at each time step is equal to the the constant drift plus the random variation of the Wiener process times the volatility term.

Once estimated the proper values for \(\mu\) and \(\sigma\) we will be ready to simulate the GBM assuming:

  • \(dt = t - (t-1) = 1\), unitary time step;
  • \(dW_t = W_t - W_{t-1} \sim \mathcal{N}(0,\ 1)\), random increments drawn from a normal distribution.

Entropy

Since its first usage in physics and thermodynamics Entropy has always been a measure of disorder and randomness of a system.

Entropy of Random Variable

Shannon expanded it defining the entropy of a random variable as the average level of information and surprise inherent to the variable’s possible outcomes through the well known equation:

\[ H(X) = \mathbb{E}[-log_b(X)] = -\sum_{x \in \mathcal{X}} p(x)\ log_b\big(p(x)\big) \]

The usual values for \(b\) are \(b=2\), \(b=e\) and, less likely, \(b=10\).

A generalization of this kind of entropy is provided by the Rényi Entropy which is defined in the following way:

\[ H_\alpha(X) = \begin{cases} \frac{1}{1-\alpha} log_b \sum_{x \in \mathcal{X}}p(x)^\alpha & \alpha > 0,\ \alpha \neq 1 \\ H(X) & \alpha = 1 \\ \end{cases} \]

Entropy of Dynamical System

Similarly to what has been done for random variables, a measure of unpredictability can be derived also for dynamical system.

Let \(x(t)\) be the trajectory of a dynamical system in the \(D\)-dimensional phase space. Divide the phase space into hypercubes of volume \(\epsilon^D\). Let \(p_{i_0, \dots, i_n}\) be the probability that the trajectory is in hypercube \(i_j\) at time \(t=j\tau\) with \(j=0,\dots, n\), where \(\tau\) is the sampling period:

\[ p_{i_0, \dots, i_n} = \mathbb{P} \big\{ x(t = j\tau) \in i_j \ \forall \ j =0, \dots, n \big\} \]

Then we define

\[ K_n = - \sum_{i_0, \dots, i_n} p_{i_0, \dots, i_n} log \ p_{i_0, \dots, i_n} \] as the amount of information needed to locate the system on a given trajectory and \(K_{N+1} - K_N\) as the information needed to predict which hypercube the trajectory will be in at time \((n+1) \tau\) given trajectories up to \(n \tau\).

The Kolmogorov entropy is then defined by:

\[ \begin{align} K & \equiv \lim_{\tau \rightarrow 0} \lim_{\epsilon \rightarrow 0} \lim_{N \rightarrow \infty} \frac{1}{N\tau} \sum_{n=0}^{N-1} \big( K_{n+1} - K_n \big) \\ & = - \lim_{\tau \rightarrow 0} \lim_{\epsilon \rightarrow 0} \lim_{N \rightarrow \infty} \frac{1}{N\tau} \sum_{i_1, \dots, i_N} p_{i_0, \dots, i_n} log \ p_{i_0, \dots, i_n} \end{align} \]

It is a crucial quantity for the characterization of chaotic systems, for instance:

  • \(K = 0\) in regular systems;
  • \(K > 0\) in chaotic systems;
  • \(K = +\infty\) (unbounded) in random system.

The Kolmogorov entropy is also related to the positive Lyapunov exponents \(\lambda_+\) of the system:

\[ K = \int \sum_+ \rho(X)\lambda_+(x)dX \]

where \(\rho(X)\) denotes the density of the attractor.

The Estimation of the Kolmogorov Entropy

The Kolmogorov Entropy offers information on how chaotic, or possibly random, the underlying system is, however, it is very difficult to obtain estimates of \(K\) without the knowledge of the system’s differential equations, directly from a time signal. We can then define a new measure of entropy based on Rényi’s entropy of order \(\alpha\):

\[ K^{(\alpha)} = - \lim_{\tau \rightarrow 0} \lim_{\epsilon \rightarrow 0} \lim_{N \rightarrow \infty} \frac{1}{N\tau} \frac{1}{\alpha - 1} log \sum_{i_1, \dots, i_N} p_{i_0, \dots, i_n}^\alpha \] It is easy to see that for \(\alpha \rightarrow 1\) it converges to the original Kolmogorov entropy \(K^{(1)} = K\) and \(K^{(q)} <= K^{(s)}\) for every \(s > q\).

The case \(\alpha=2\), \(K^{(2)}=K_2\) from now on, is of particular interest because of the following properties:

  1. \(K_2 \geq 0\)
  2. \(K_2 \leq K\), it is a lower bound on the true Kolmogorov entropy
  3. \(K_2 \neq 0\) for chaotic systems
  4. \(K_2 \rightarrow \infty\) for random systems

and it turns out to be also numerically close to \(K\) in many typical cases.

However, the most important property of \(K_2\) is that it can be easily estimated for experimental signals using the generalized correlation integrals \(C_d(r)\):

\[ \begin{align} C_{d}(r) &= \lim_{N \rightarrow \infty} \frac{1}{N^2} \bigg[ \text{Numer of pairs } i, j \text{ with } \Big( \sum_{k=1}^{d} |X_{i+k} - X_{j+k}|^2 \Big)^{1/2} < r \bigg] \\ &= \lim_{N \rightarrow \infty} \frac{1}{N^2} \sum_{i=1}^N \sum_{j=1}^N \Bigg( H \bigg( r - \Big( \sum_{k=1}^{d} |X_{i+k} - X_{j+k}|^2 \Big)^{1/2} \bigg) \Bigg) \end{align} \] where \(H\) is the Heaviside function, \(r\) is the distance threshold and \(d\) is the number of measurements used for each coordinate.

It can be shown that \(C_d(r) \underset{r \rightarrow 0 \\ d \rightarrow \infty}{\sim} r^\nu exp(-d \tau K_2)\) where \(\nu\) denotes the correlation exponent, moreover, it has been proved that \(\nu \lesssim D\) where \(D\) is the fractal dimension of the attractor.

From the last expression we can finally derive

\[ \begin{align} &K_{2,d}(r) = \frac{1}{\tau} log \frac{C_d(r)}{C_{d+1}(r)} & \lim_{r \rightarrow 0 \\ d \rightarrow \infty} K_{2,d}(r) \sim K_2 \end{align} \]

Hurst Exponent

The Hurst exponent is a measure widely used to characterize time series, it quantifies the relative tendency of a time series either to regress strongly to the mean or to cluster in a direction.

The values of the Hurst exponent range between \(0\) and \(1\) and based on the value of \(Hu\), we can classify any time series into one of the three categories:

It is defined in terms of the asymptotic behavior of the rescaled range:

\[ \mathbb{E} \bigg[ \frac{R(n)}{S(n)} \bigg] = Cn^{Hu} \text{ as } n \rightarrow \infty \] where:

Simulation

In this section we will simulate two types of GBM: with drift (\(\mu=\bar{\mu}\) where \(\bar{\mu}\) is the SP500 average daily return) and without drift (\(\mu=0\)). In both cases the volatility \(\sigma=\bar{\sigma}\) will be set equal to the SP500 average daily volatility.

In order to do that first we have to compute the average return \(\bar{\mu}\) and the average volatility \(\bar{\sigma}\).

Let \(S_t\) be the stock price, we define \(PC(t) = \frac{S_{t+1} - S_t}{S_t}\) the percentage change between two time instants and \(R(t) = log \ S_{t+1} - log \ S_t = log(1 + PC(t))\) as the log-return.

We can now proceed with the estimation of \(\bar{\mu}\) and \(\bar{\sigma}\), respectively:

\[ \begin{align} & \bar{\mu} = \mathbb{E} \big[ PC(t) \big] & \bar{\sigma}^2 = \mathbb{V}ar \big[ PC(t) \big] \end{align} \]

The data used are the daily SP500 prices from 1950-01-03 to 2015-12-31.

# Evaluate mean and standard deviation of percentage returns
sp500_price <- rev(sp500$open)
pct_change_mean <- mean(diff(sp500_price) / sp500_price[-length(sp500_price)])
pct_change_sd <- sd(diff(sp500_price) / sp500_price[-length(sp500_price)])
print(paste('Mean Pct-Returns', format(pct_change_mean, scientific = T), 
            ', Sd Pct-Returns', format(pct_change_sd, scientific = T)))
## [1] "Mean Pct-Returns 3.355782e-04 , Sd Pct-Returns 9.505787e-03"
# Define mu_bar and sigma_bar
mu_bar <- pct_change_mean
sigma_bar <- pct_change_sd

To simulate a GBM of length \(n\) we have to:

  1. evaluate \(n-1\) percentage changes \(PC(t) = \mu + \sigma w(t) \ \forall \ t=1, \dots, (n-1)\), where \(w(t)\) are iid \(\mathcal{N}(0, 1)\);
  2. compute GBM signal according to \(S(t) = S_0 \prod_{i=1}^t PC(t)\).
simulate_GBM <- function(n, mu, sigma, S0 = 1){
  w <- rnorm(n-1, 0, 1)
  pc <- mu + sigma * w
  signal <- cumprod(c(S0, 1+pc))
  return(signal)
}

GBM Without Drift \(\mu = 0\)

First we analyze the case without drift, \(\mu = 0\).

n <- 5000
mu <- 0
sigma <- sigma_bar
gbm <- simulate_GBM(n, mu, sigma)

Here we simulate several GBM and we plot them on both linear and log scale.

Kolmogorov Entropy Estimation - \(K_2\)

Unfortunately there is no library in R implementing the \(K_2\) entropy therefore we had to rely on the python package EntropyHub to perform such estimates.

Here we compare the \(K_2\) entropy estimated for several signals, specifically:

  • the GBM simulated above \(S(t)\)
  • the log-returns for the GBM simulated above \(R(t)\)
  • a periodic signal \(Sine(t)\)
  • a noise signal \(Noise(t)\) (normal iid samples)
import EntropyHub
import numpy as np
from matplotlib import pyplot as plt
import pandas as pd
import seaborn as sns

_s = np.array(r.gbm)
_r = np.diff(np.log(_s))
_sine = np.sin(np.linspace(0, 10, _s.size)) * np.std(_s)
_noise = np.random.normal(size=_s.size) * np.std(_s)

fig, ((ax1, ax2), (ax3, ax4)) = plt.subplots(nrows=2, ncols=2)

ax1.plot(_s, label='S(t)', alpha=.8)
ax2.plot(_r, label='R(t)', alpha=.8)
ax3.plot(_sine, label='Sine(t)', alpha=.8)
ax4.plot(_noise, label='Noise(t)', alpha=.8)
ax1.legend()
ax2.legend()
ax3.legend()
ax4.legend()
plt.show()

d = 15
_s_entropy = EntropyHub.K2En(_s, m=d)[0]
print('S(t) entropy:', _s_entropy)
## S(t) entropy: [0.3753156  0.28019882 0.2832476  0.30459839 0.32618082 0.34475932
##  0.35956604 0.3761557  0.38318545 0.38832004 0.39215748 0.39119164
##  0.38812988 0.38183007 0.36723231]
_r_entropy = EntropyHub.K2En(_r, m=d)[0]
print('R(t) entropy:', _r_entropy)
## R(t) entropy: [2.43744498 2.60022125 2.71790524 2.80464207 2.48450617        nan
##         nan        nan        nan        nan        nan        nan
##         nan        nan        nan]

Increasing the value of \(d\) makes the procedure numerically unstable, for this reason we will use \(d=5\).

d = 5
_s_entropy = EntropyHub.K2En(_s, m=d)[0]
_r_entropy = EntropyHub.K2En(_r, m=d)[0]
_sine_entropy = EntropyHub.K2En(_sine, m=d)[0]
_noise_entropy = EntropyHub.K2En(_noise, m=d)[0]

From the plot above some considerations can be carried out:

  • evaluate the entropy using \(R(t)\) instead of \(S(t)\) makes the results more reliable;
  • due to numerical issues it is not possible to discriminate between chaotic and random processes, however;
  • the comparison with a noise process can give us insight on how a random process would behave.

Hurst Exponent

Here we test the hurst exponent of a GBM without drift.

M = 50
# Hurst exponent of S(t)
hu_s_df <- data.frame(t(sapply(1:M, function(x) unlist(hurstexp(simulate_GBM(n, mu, sigma), display = F)))))
melt(hu_s_df) %>%
  ggplot() + 
  geom_histogram(aes(x=value), bins = 15) + 
  ggtitle('Hurst Exponent of S(t)') +
  facet_wrap( ~ variable)

# Hurst exponent of R(t)
hu_r_df <- data.frame(t(sapply(1:M, function(x) unlist(hurstexp(diff(log(simulate_GBM(n, mu, sigma))), display = F)))))
melt(hu_r_df) %>%
  ggplot() + 
  geom_histogram(aes(x=value), bins = 15) + 
  ggtitle('Hurst Exponent of R(t)') +
  facet_wrap( ~ variable)

The Hurst exponent shows something we already noted before analyzing the \(K_2\) values: when computed on the GBM the results are unreliable. For instance, when computed on the GBM itself almost all the \(Hu\) estimates agree on the signal being strongly trending while, when computed using the log-returns, they show uncertainty on the trending/mean reverting property of the signal.

GBM With Drift \(\mu =\) 3.355782e-04

Here we repeat the same steps we performed before and we analyze if the presence of a drift component produce any difference.

n <- 5000
mu <- mu_bar
sigma <- sigma_bar
gbm <- simulate_GBM(n, mu, sigma)

Here we simulate several GBM and we plot them on both linear and log scale.

Kolmogorov Entropy Estimation - \(K_2\)

As before we compare the \(K_2\) entropy estimated for several signals, specifically:

  • the GBM simulated above \(S(t)\)
  • the log-returns for the GBM simulated above \(R(t)\)
  • a periodic signal \(Sine(t)\)
  • a noise signal \(Noise(t)\) (normal iid samples)
_s = np.array(r.gbm)
_r = np.diff(np.log(_s))
_sine = np.sin(np.linspace(0, 10, _s.size)) * np.std(_s)
_noise = np.random.normal(size=_s.size) * np.std(_s)

fig, ((ax1, ax2), (ax3, ax4)) = plt.subplots(nrows=2, ncols=2)

ax1.plot(_s, label='S(t)', alpha=.8)
ax2.plot(_r, label='R(t)', alpha=.8)
ax3.plot(_sine, label='Sine(t)', alpha=.8)
ax4.plot(_noise, label='Noise(t)', alpha=.8)
ax1.legend()
ax2.legend()
ax3.legend()
ax4.legend()
plt.show()

d = 5
_s_entropy = EntropyHub.K2En(_s, m=d)[0]
_r_entropy = EntropyHub.K2En(_r, m=d)[0]
_sine_entropy = EntropyHub.K2En(_sine, m=d)[0]
_noise_entropy = EntropyHub.K2En(_noise, m=d)[0]

Looking a the plots above we can formulate the following considerations:

  • as before, evaluate the entropy using \(R(t)\), instead of \(S(t)\), makes the results more reliable;
  • the GBM signal \(S(t)\) shows an higher predictability, while;
  • the log-returns signal \(R(t)\) is indistinguishable from the random process.

Hurst Exponent

Now we test the Hurst exponent of a GBM with drift. In this case we expect that the drift acts as the trend and therefore we’re expecting that, on average, \(Hu > 0.5\).

M = 50
# Hurst exponent of S(t)
hu_s_df <- data.frame(t(sapply(1:M, function(x) unlist(hurstexp(simulate_GBM(n, mu, sigma), display = F)))))
melt(hu_s_df) %>%
  ggplot() + 
  geom_histogram(aes(x=value), bins = 15) + 
  ggtitle('Hurst Exponent of S(t)') +
  facet_wrap( ~ variable)

# Hurst exponent of R(t)
hu_r_df <- data.frame(t(sapply(1:M, function(x) unlist(hurstexp(diff(log(simulate_GBM(n, mu, sigma))), display = F)))))
melt(hu_r_df) %>%
  ggplot() + 
  geom_histogram(aes(x=value), bins = 15) + 
  ggtitle('Hurst Exponent of R(t)') +
  facet_wrap( ~ variable)

Despite what we expected, the \(Hu\) estimates show uncertainty on the trending/mean reverting property of the signal. We think that a greater \(\mu / \sigma\) ratio would help making \(Hu > 0.5\): a stronger trend in conjunction with a smaller volatility.

The Stock Market

Now we want to use what we have defined until this point to investigate the stock market.

The Data

The data we are going to use are the daily closing price of \(441\) stocks going from 2010-01-05 to 2021-12-27 (3018 samples for each stock).

data <- read.csv('prices.csv')
knitr::kable(head(data[, 1:6]))
Date A AAL AAP AAPL ABC
2010-01-05 20.23958 5.005957 38.27243 6.564354 22.01566
2010-01-06 20.16766 4.798553 38.60616 6.459941 21.80749
2010-01-07 20.14152 4.939965 38.59660 6.447997 21.45777
2010-01-08 20.13498 4.845692 38.74916 6.490865 21.69092
2010-01-11 20.14806 4.751417 38.36779 6.433607 21.93239
2010-01-12 19.90617 4.789125 37.70034 6.360424 22.08227

Each price time series is divided by the price in the first time instant.

prices <- apply(data[, -1], 2, function(x) x/x[1])

Kolmogorov Entropy Estimation - \(K_2\)

We can now select \(M=50\) stocks at random and compute their \(K_2\) values using the price signal \(S(t)\). We show both the raw and the aggregated \(K_2\) values comparing the latter with the values we got before for a noise and a cyclical signal.

import seaborn as sns
_prices = r.prices
d = 5
ncol = _prices.shape[1]
M = 50
_assets = np.random.choice(ncol, M)
K2_prices = np.zeros((M, d))
for i, idx in enumerate(_assets) :
  K2_prices[i] = EntropyHub.K2En(_prices[:, idx], m=d)[0]
  
K2_prices_df = pd.DataFrame(K2_prices, columns = range(1, d+1)).reset_index().melt('index')
K2_prices_df.columns = ['Stock Index', 'd', 'K2']

Here we do the same using the log-returns \(R(t)\) instead.

_log_returns = np.diff(np.log(_prices), axis=0)

K2_log_returns = np.zeros((M, d))
for i, idx in enumerate(_assets) :
  K2_log_returns[i] = EntropyHub.K2En(_log_returns[:, idx], m=d)[0]
  
K2_log_returns_df = pd.DataFrame(K2_log_returns, columns = range(1, d+1)).reset_index().melt('index')
K2_log_returns_df.columns = ['Stock Index', 'd', 'K2']

The plots strongly resemble the ones we got before, this suggests that the prices may indeed follow a GBM with drift.

Hurst Exponent

Now we proceed evaluating the Hurst Exponent \(Hu\) for all the assets’ log-returns. We have indeed already shown that the evaluate \(Hu\) directly from the price signal often leads to incorrect results.

price_log_returns <- diff(log(prices))
# Hurst exponent of R(t)
hu_price_log_returns_df <- data.frame(t(apply(price_log_returns, 2, function(x) unlist(hurstexp(x, display = F)))))
melt(hu_price_log_returns_df) %>%
  ggplot() + 
  geom_histogram(aes(x=value), bins = 15) + 
  ggtitle('Hurst Exponent of R(t)') +
  facet_wrap( ~ variable)

From the plot above we can clearly see that many of the \(Hu\) distributions are centered in \(Hu = 0.5\), value such that the process is considered a random walk.

Conclusion

In this project we have investigated the stock market through the usage of the Kolmogorov Entropy and the Hurst Exponent.

Though the assumption of the price being a Geometric Random Walk can be proven false, i.e. looking at the Autocorrelation Function of any random stock below,

set.seed(0)
stock <- sample(colnames(prices),1)
print(paste('Stock:', stock))
## [1] "Stock: UDR"

our experiments gave us many insights on the unpredictable nature of the stock market. Often prices and random processes were indistinguishable according to the measures we presented, see the distributions of the Hurst exponents \(Hu\) centered around \(0.5\) or the \(K_2\) values close to the noise level.

We believe that the unpredictability of this system is rooted in amount of interactions it has to deal with, such complexity makes the prediction of the following price impossible. Studies in the field of high frequency trading have shown that the prediction window, time span within which is possible to perform prediction with high accuracy, lasts, indeed, just few millisecond.